# R 包与数据导入
library(Seurat)
library(circlize)
library(tidyverse)
library(SCENIC)
library(RcisTarget)
library(AUCell)
library(gplots)
library(ComplexHeatmap)
library(patchwork)
library(RColorBrewer)
load("demo_data/project_scenic.rdata")

# 初始化 SCENIC
org <- "hgnc"
dbDir <- "demo_data/database/human/"
dbs <- defaultDbNames[[org]]
data(list="motifAnnotations_hgnc_v9", package="RcisTarget")
motifAnnotations_hgnc <- motifAnnotations_hgnc_v9
scenicOptions <- initializeScenic(org = org, dbDir = dbDir, dbs = dbs)

# 提取原始表达矩阵，并执行过滤
exprMat <- as.matrix(GetAssayData(project,slot='counts',assay='RNA'))
genesKept <- geneFiltering(exprMat,
                           scenicOptions = scenicOptions,
                           minCountsPerGene = 6,
                           minSamples = ncol(exprMat)*0.01)
exprMat_filtered <- exprMat[genesKept,] 

# 构建共表达网络
scenicOptions@settings$seed <- 100              # 设置随机数种子，保证每次运行的结果都是一致的。
scenicOptions@settings$nCores <- 10             # 设置线程数，根据服务器性能进行设置
# 注意：不要在个人电脑上运行 GENIE3
# 方案一
runGenie3(exprMat_filtered, scenicOptions)      # 运行 GENIE3 构建共表达网络

# 替代方案
# 方案二
# 使用 pyscenic 的 pyscenic grn 来代替 runGenie3，仅仅是这步使用 pyscenic 进行替代，其余步骤与 R 版本的 scenic 一致。
data2grn <- exprMat_filtered %>% as.data.frame() %>% tibble::rownames_to_column(var = "Symbol")
data.table::fwrite(data2grn, file = 'matrix_for_pyscenic.csv', quote = FALSE, row.names = FALSE, sep = ",")
system("pyscenic grn --method grnboost2 --output ./step1_mx_adjacencies.tsv --num_workers 20 --seed 100 --transpose ./matrix_for_pyscenic.csv /public/home/lxw/05.Test/SCENIC/database/human/hs_hgnc_tfs.txt")
grn_result <- data.table::fread("step1_mx_adjacencies.tsv") %>% as.data.frame()
#     TF target importance
#1  SPI1 TYROBP   133.9828
#2  RPL6  RPL34   115.3732
#3  RPL6   RPS6   108.9242
# GENIE3的结果文件列名与 pyscenic grn运行结果的文件列名不一致，改一下，保存到 int/1.4_GENIE3_linkList.Rds，接下来的步骤不变。
# 经测试 num_workers 为 20 的情况下半小时内就能构建完共表达网络，相较于 runGenie3 速度快了不止一星半点
colnames(grn_result) <- c("TF","Target","weight")
saveRDS(grn_result,file = "int/1.4_GENIE3_linkList.Rds")

runCorrelation(exprMat_filtered, scenicOptions)
# 构建基因调控网络并进行打分(Build and score the GRN)，运行过程很耗时，不要在个人电脑上尝试
scenicOptions <- runSCENIC_1_coexNetwork2modules(scenicOptions)
scenicOptions <- runSCENIC_2_createRegulons(scenicOptions, coexMethod=c("top10perTarget")) 
scenicOptions <- runSCENIC_3_scoreCells(scenicOptions, log2(exprMat_filtered + 1))

# 读入已分析完成的结果
regulonAUC <- readRDS("results/int/3.4_regulonAUC.Rds")
# onlyNonDuplicatedExtended 它是 SCENIC 自带的函数，当某一个 regulon 既存在高质量，又高质量时，对低质量(带_extended后缀)进行过滤；
# 当某一个 regulon 只有低质量时，不进行过滤。
regulonAUC <- regulonAUC[onlyNonDuplicatedExtended(rownames(regulonAUC)),]
regAct <- getAUC(regulonAUC)
# 带有 _extended 后缀的都属于 低置信度 的 regulon，在后续的分析应过滤
regAct <- regAct[!str_detect(rownames(regAct),"extended"),]

# 单细胞水平热图展示
regAct <- regAct[apply(regAct, 1, sd) != 0, ]
regAct <- t(scale(t(regAct), center = T, scale=T))
column_df <- data.frame(CellType=project@meta.data$CellType, Sample=project@meta.data$Sample)
celltype_color <- colorRampPalette(brewer.pal(9, 'Set1'))(length(levels(project@meta.data$CellType)))
names(celltype_color) <- levels(project@meta.data$CellType)
sample_color <- colorRampPalette(brewer.pal(9, 'Set1'))(length(levels(project@meta.data$Sample)))
names(sample_color) <- levels(project@meta.data$Sample)
col_fun <- circlize::colorRamp2(seq(-5,5,0.2),gplots::colorpanel(51, low="blue", mid="white",high="red"))
column_annotation <- HeatmapAnnotation(df = column_df, col = list(Cluster = celltype_color, Sample = sample_color))
# 绘图
complexheatmap <- Heatmap(regAct, 
                          col = col_fun, 
                          name = "heatmap", 
                          heatmap_legend_param = list(legend_direction = "horizontal", 
                                                      legend_width = unit(10, "cm"), 
                                                      title_position = "lefttop"),
                          column_names_side = "top", 
                          cluster_columns=TRUE, 
                          cluster_rows=TRUE, 
                          show_column_names = FALSE, 
                          row_names_gp = gpar(fontsize = 12), 
                          top_annotation = column_annotation)
# 保存图片
png(file="SCENIC_AUC_score_heatmap.png", width = 20, height = 15, units = "in",res = 300)
draw(complexheatmap, heatmap_legend_side = "bottom")
dev.off()


# 细胞亚群找差异(特异)的 TF 及 热图展示
regulonAUC <- readRDS("int/3.4_regulonAUC.Rds")
regAct <- getAUC(regulonAUC)
regAct <- regAct[!str_detect(rownames(regAct),"extended"),]
regAct_t <- t(regAct)
regAct_t <- data.frame(CellType = project@meta.data$CellType, regAct_t)
gene_sets <- colnames(regAct_t)[-1]
celltype <- unique(project@meta.data$CellType)
score_tscore <- matrix(data=NA, nrow=length(gene_sets), ncol=length(celltype))
score_pvalue <- matrix(data=NA, nrow=length(gene_sets), ncol=length(celltype))
rownames(score_tscore) <- rownames(regAct)
colnames(score_tscore) <- celltype
rownames(score_pvalue) <- rownames(regAct)
colnames(score_pvalue) <- celltype

for(index in 1:length(celltype)){
  tv <- NULL
  pv <- NULL
  for(a_gene_set in gene_sets){
    ttest_re <- t.test(regAct_t[, a_gene_set][regAct_t$CellType == celltype[index]], 
                       regAct_t[, a_gene_set][regAct_t$CellType != celltype[index]])
    tv <- c(tv, ttest_re$statistic[[1]])
    pv <- c(pv, ttest_re$p.value)
  }
  score_tscore[, index] <- tv
  score_pvalue[, index] <- pv
}

column_df <- data.frame(CellType = celltype)
celltype_color <- colorRampPalette(brewer.pal(9, 'Set1'))(length(celltype))
names(celltype_color) <- celltype
column_annotation <- HeatmapAnnotation(df = column_df, col = list(CellType = celltype_color))
# 绘图
complexheatmap <- Heatmap(score_tscore, 
                          col=colorpanel(75, low="blue", mid="white",high="red"), 
                          name = "heatmap", 
                          heatmap_legend_param = list(legend_direction = "horizontal", 
                                                      legend_width = unit(10, "cm"), 
                                                      title_position = "lefttop"),
                          column_names_side = "top", 
                          cluster_columns=TRUE, 
                          cluster_rows=TRUE, 
                          show_column_names = FALSE, 
                          row_names_gp = gpar(fontsize = 12), 
                          top_annotation = column_annotation)
# 保存图片
png(file="SCENIC_AUC_score_celltype_heatmap.png", width = 20, height = 15, units = "in",res = 300)
draw(complexheatmap, heatmap_legend_side = "bottom")
dev.off()

# regulons 调控得分的降维图
regulonAUC <- readRDS("int/3.4_regulonAUC.Rds")
regulonAUC <- regulonAUC[onlyNonDuplicatedExtended(rownames(regulonAUC)),]
regAct <- getAUC(regulonAUC)
regAct <- regAct[!str_detect(rownames(regAct),"extended"),]
project[["SCENIC"]] <- CreateAssayObject(counts = regAct)
DefaultAssay(project) <- "SCENIC"
regulons_list <- c("RUNX3 (168g)")
p1 <- FeaturePlot(object = project, features = regulons_list, reduction='umap', pt.size=0.1)
p2 <- FeaturePlot(object = project, features = regulons_list, reduction='tsne', pt.size=0.1)
p <- p1 | p2
# 保存图片
ggsave(p, filename = "tsne_umap_SUNX10.png",width = 10,height = 5,dpi = 300)

# regulons 调控山脊图和小提琴图
p3 <- RidgePlot(project, features = regulons_list, group.by = "CellType") + 
  labs(x = NULL,y = NULL) + 
  theme(axis.line = element_line(colour="black"),
        panel.grid = element_blank(),
        legend.position = 'none',
        plot.margin = margin(c(20,20,20,20),"mm"),
        strip.background = element_rect(colour="#f0f0f0",fill="#f0f0f0"),
        strip.text = element_text(face="bold"))

p4 <- VlnPlot(project, features = regulons_list, group.by = "CellType", pt.size = 0) + 
  labs(x = NULL,y = "AUC score") + 
  theme(axis.line = element_line(colour="black"),
        panel.grid = element_blank(),
        legend.position = 'none',
        plot.margin = margin(c(20,20,20,20),"mm"),
        strip.background = element_rect(colour="#f0f0f0",fill="#f0f0f0"),
        strip.text = element_text(face="bold"))
p <- p3 | p4
# 保存图片
ggsave(p, filename = "vlnplot_ridgeplot_SUNX10.png",width = 10,height = 5,dpi = 300)